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We study the isotope effect on the temperature of the proton order/disorder phase transition 
between ice XI and ice Ih, using the quasiharmonic approximation combined with ab initio density 
functional theory calculations. We show that this method is accurate enough to obtain a phase 
transition temperature difference between light ice (H 2 O) and heavy ice (D 2 O) of 6 K as compared 
to the experimental value of 4 K. More importantly, we are able to explain the origin of the isotope 
effect on the much debated large temperature difference observed in the phase transition. The source 
of the difference is directly linked to the physics behind the anomalous isotope effect on the volume 
of hexagonal ice that was recently explained in [Phys. Rev. Lett. 108, 193003 (2012)]. These 
results indicate that the same physics might be behind the isotope effects in transition temperatures 
between other ice phases. 


I. INTRODUCTION 

The polymorphism of ice is revealed by its rich phase 
diagramiii^. The availability of different proton configu¬ 
rations that satisfy Bernal-Fowler “ice-rules’— adds an¬ 
other dimension to this phase diagram, given that the 
same crystalline structure could exist in proton-ordered 
and disordered form. This leads to additional phases, 
separated by their corresponding order/disorder phase 
transitions, as in the case of ice Xl/ice Ih, ice IX/ice III, 
and ice Vlll/ice VII^. 

In this paper we focus on the phase transition be¬ 
tween proton-ordered (ice XI) and proton-disordered (ice 
Ih) hexagonal ice. This phase transition has been sub¬ 
ject of a large number of experimental^^— and theoretical 
studies^Ti^i. However, open questions remain about the 
mechanisms behind the phase transition and the impor¬ 
tance of nuclear quantum effects in the temperature of 
the transitions^. 

Experimentally, it is difficult to observe the phase tran¬ 
sition from ice Ih to ice XL A glass transition occurs 
at around 100-110 KSiS, diminishing the proton mobility 
and locking protons in their disordered positions, before 
they orient to form the proton-ordered ice XI structure. 
This is overcome by catalyzing ice Ih by KOHSt^, which 
allows the lattice parameters of both proton-ordered ice 
^j23-_28 disordered ice I h^s, 29,30 experimen¬ 

tally measured. The order/disorder phase transition is 
achieved at 72 K for light^ and 76 K for heavy ice^. Al¬ 
though the isotope effect on the phase transition temper¬ 
ature is measured to be 4 K, a theoretical explanation for 
this difference is still missing. Ref. |3l| and IH associated 
the origin of the isotope effect on the transition temper¬ 
ature to the difference in vibrational energies between 
the two phases estimating a ^ 27 K isotope effect on 
the transition temperature, while Ref. sought the ex¬ 
planation in the difference in reorientation of the dipole 
moments of heavy and light ices and predicted a much 


smaller ~ 1 K isotope effect. However, in both of these 
studies they assume no isotope effect on the volume of 
the two ices and they also did not take into account the 
competing anharmonicities between the intra-molecular 
covalent bonds and the inter-molecular hydrogen bonds. 
In this work, we reexamine the issue taking into account 
these two additional effects, extending a previous work 
on the anomalies in the isotope effect on the volume of 
ice2^. 

A large literature is devoted to the study of the or¬ 
der/disorder phase transition in hexagonal ice. Two main 
questions are discussed, (i) the ordering nature in the low 
temperature phase, and (ii) theory and simulation predic¬ 
tions for the phase transition temperature. Experiments 
such as neutron diffractioi>2^“— , as well as measurements 
performed under an electric field^^ indicate that the or¬ 
dered phase has ferroelectric order. However, among the¬ 
ory and simulation works there is a large dispersion of 
results and lack of agreement^r^. The predicted low T 
stable phase depends strongly on the choice of bound¬ 
ary conditions, electrostatic multipoles, and treatment of 
long range interactions— In addition, semi-empirical 
force field models fitted to reproduce experimental data 
are not accurate enough to distinguish small energy dif¬ 
ferences between different proton orderings. 

According to Ref. [^, the TIP4P-FQ24 model pre¬ 
dicts the proton-disordered phase to be stable, in agree¬ 
ment with Ref. [ 1 ^ where other less known models were 
studied. On the other hand popular water models like 
SPC/E2^, TIP4P^, TIP5P-ES and NvdBH models pre¬ 
dict the proton ordered phases to be more stable in the 
low temperature limit. Among these studies, only the 
NvdEi^ model predicts the ferroelectric-ordered phase as 
the lowest energy phase, in agreement with experiments, 
while the other three models predict the stable phase to 
be antiferroelectric-ordered. In this same study, it was 
also shown that modifying the polarizability of the KW- 
pol model was enough to favor ferroelectric ordering over 
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disordered configurations at low T— . Therefore polar¬ 
izability is an important factor in obtaining the correct 
potential energy surface. 

The effect of proton disorder on the hexagonal ice 
structure has also been studied using ab initio density 
functional theory (DFT). DFT calculations correctly re¬ 
produce the lattice structure^^, and give the cohesive en¬ 
ergy of ferroelectric-ordered ice to be larger than either 
antiferroelectric-ordered or disordered ice a^^d^ . That is, 
ferroelectric-ordered ice is the stable ground state. A 
recent DFT study of ice slabs shows that contrary to 
the bulk case, the antiferroelectric-ordered ice is more 
stable in the case of thin films^. DFT-based simula¬ 
tions have also been used to predict the phase transition 
temperature. A Monte Carlo study, where DFT calcu¬ 
lations of hydrogen bond configuration energies are used 
to parametrize a model to perform Monte Carlo simula¬ 
tions, predicted ice XI as the most stable phase, with 
a transition temperature of 98 K.— Another recent 
DFT-based Monte Carlo study of dielectric properties 
of ice, predicted the temperature of the order-disorder 
phase transition to be around 70-80 K— . The advan¬ 
tage of DFT-based Monte Carlo simulations is that they 
can explore the configurational entropy of the free en¬ 
ergy surface in good detail. However, none of these cal¬ 
culations include zero point nuclear quantum effects, or 
thereby investigate the transition temperature difference 
between different isotopes. Our goal is to investigate the 
order/disorder transition from ferroelectric-ordered ice or 
antiferroelectric-ordered ice to disordered ice, with differ¬ 
ent isotopic compositions, from an ab initio perspective, 
including nuclear quantum effects. 

In our study we will also compare a polarizable force 
field model, TTM3-F— to DFT calculations for the pre¬ 
diction of the most stable phase at low temperatures in¬ 
cluding zero point corrections. 


II. THEORY 


In a recent study, we explained the anomalous isotope 
effect on the volume of ic o^^i^°'^^ , by obtaining the free 
energy with ab initio DFT within the quasiharmonic ap¬ 
proximation. We have shown that the anticorrelation 
between the intra-molecular OH covalent bonds and the 
inter-molecular hydrogen bonds makes the volume per 
molecule of D 2 O ice larger than that of H 2 O ice^. 

In this work, we extend our study of nuclear quan¬ 
tum effects to analyze the contribution to the or¬ 
der/disorder phase transition using both ab initio 
DFT functionals and the TTM3-F— force field model. 
We investigate both ferroelectric-ordered/disordered and 
antiferroelectric-ordered/disordered phase transitions. 
In addition, we analyze the importance of van der Waals 
forces, by comparing a generalized gradient approxi¬ 
mated functional, Perdew-Burke-Ernzerhof (PBE)iS, to 
a van der Waals functiona l"^ vdW-DF^®® . We ob¬ 
tain the temperature dependence of the free energy for 


both ice phases using the quasiharmonic approximation 
(QHA), and we compare the phase transition tempera¬ 
ture of different isotopes. 


A. Free Energy within Quasiharmonic 
Approximation 

To account for nuclear quantum effects, quantum har¬ 
monic eigenstates are needed as a function of volume V, 
at volumes near Vq, the “frozen lattice” zero pressure 
volume that minimizes the Born-Oppenheimer energy, 
Eo{V). To lowest order in a Taylor series around Vq, 
we have 

EoiV)=EoiVo) + ^^{V-Vor ( 1 ) 

and 

ujk{V) = oj{Vo) (^1 - (2) 

Bq is the dominant part of the bulk modulus, omitting 
vibrational corrections which will be discussed in a later 
paper. The “mode Griineisen parameters” 'jk are defined 
as 


cl(lna;fe) _ V dujk 
9(lnE) “ oJkdV 


( 3 ) 


The phonon frequencies, ujk are calculated at three differ¬ 
ent volumes. The volume dependence of uJk(V) is calcu¬ 
lated to the linear order. Then the Helmholtz free energy 
of independent harmonic oscillators acquires a 
volume-dependence through uJk{V), 


FiV,T) = EoiV) + 

'hiokiv) 


E 


k ^ 

-TSh 


+ fc^rin (1 - e-^‘^k{v)/kBT^^ 


( 4 ) 


The index k runs over both phonon branches and phonon 
wave vectors within the Brillouin zone. This “quasihar¬ 
monic” approximation is correct to first order for volume 
derivatives like P = —{dF/dV)T- Higher volume deriva¬ 
tives, such as B{T), in general may require higher vol¬ 
ume derivatives of Eq and uJk- As shown in our recent 
contributions^^iii, the first derivative in eq. [2] is a good 
approximation for hexagonal ice. The temperature de¬ 
pendence of volume Vp^i^iT) is then found in the usual 
way by minimizing F{V,T) at fixed T, the same as set¬ 
ting P{T) = 0. 

The last part of the free energy, Sh, is the entropy 
of the proton disorder. This term is zero for proton- 
ordered ice phases. For the proton-disordered phase, ice 
Ih, we use the estimation by Pauling, Sh = Nks ln(3/2), 
which was obtained by counting hydrogen orientations 
that obey the ice rules^, and experimentally confirmed 
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for fully disordered case a^^d^ . We assume that this term 
does not change with temperature. 

Lastly, the classical limit of the free energy is obtained 
by taking the high temperature limit of the QHA: 


F{V,T) = Eo{V) + 


E 


fcsTln 


V keT ) 


-TSh. 


(5) 




B. Cohesive Energy 


To determine which structure is the most stable one 
at zero temperature, cohesive energies of ices are calcu¬ 
lated without (Ac) s-iid with (Ec) zero point effects. The 
cohesive energy is defined as the amount the energy of a 
molecule is lowered in a crystal relative to in vacuum: 


/T'lce 

^0 _ ~^0 _ _ ^monomer 

^molecules 


( 6 ) 


Ac 


Amolecules 


TTimonomer 77’monomer 

^0 ^vib 


(7) 


where the vibrational energy of the three 

modes of the monomer is A™™°™®''. The classical cohe¬ 
sive energy, E^ is defined in eq. [B] using the Kohn-Sham 
energies of the ice and monomer; and similarly, the quan¬ 
tum cohesive energy with zero point effects, Ac is defined 
in eq. [T] 


FIG. 1. Unit cell of the ferroelectric proton-ordered ice XI 
structure. The 4 molecules in the unit cell are labeled with a 
star symbol next to it, and a, and c lattice vectors are shown. 
The image on the left is the side view of the x-z plane; the 
image on the right is the top view of the x-y plane. 


moment, as shown in Fig. [5] Similarly, we have used a 
unit cell of 8 molecules for the DFT, and 96 molecules 
for the TTM3-F calculations. 


4. 4. V' 


I 







FIG. 2. Antiferroelectric proton-ordered ice aXI structure. 
The 8 molecules in the unit cell are labeled with a star symbol 
next to it, and a, b, and c lattice vectors are shown. The image 
on the left is the side view of the x-z plane; the image on the 
right is the top view of the x-y plane. 


III. SIMULATION DETAILS 
A. System Description 

In order to predict the most stable phase of hexago¬ 
nal ice in the zero temperature limit, we performed total 
energy calculations of three hexagonal ices with different 
proton configurations, (i) Ice XL This is the ferroelectric- 
proton-ordered ice. Oxygen atoms are constrained to the 
hexagonal wurtzite lattice and hydrogen atoms are or¬ 
dered such that ice XI has a net dipole moment along 
the c axis, shown in Fig. [TJ Precise measurements of lat¬ 
tice structure of ice XI have shown ferroelectric ordering, 
with a net dipole moment along the c-axisi^“— There 
are 4 molecules per formula unit. However, TTM3-F cal¬ 
culations were performed for a 3a x 2-\/3a x 2c supercell 
with 96 molecules, with the same cell size as disordered 
ice Ih. 

(ii) Ice aXI. We label the antiferroelectric-proton- 
ordered ice as ice aXI. Ice aXI has 8 molecules in the 
unit cell. The unit cell is doubled from the ferroelectric 
proton-ordered ice XI along the x-y plane, with dipole 
moments of the neighboring molecules pointing in oppo¬ 
site directions such that the system has no net dipole 
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FIG. 3. Proton-disordered ice Ih structure. The image on the 
left is the side view of the x-z plane; the image on the right 
is the top view of the x-y plane. 

(iii) Ice Ih. Experimentally, the lattice structure of 
light(H 20 ) and heavy(D 20 ) proton-disordered hexago¬ 
nal ice Ih have been measured using both syncrotron 
radiation^Si^ and neutron diffraction^^, with good agree¬ 
ment. Oxygen atoms still have an underlying hexagonal 
lattice, while hydrogen atoms are disordered such that it 
has no net dipole moment. An example of this system 
is shown in Fig. [3] To accommodate different proton- 
disordered configurations in ice Ih, we have used a large 
cell of 96 molecules with dimension 3a x 2-\/3a x 2c in our 
DFT calculations. 
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We have computed five diflFerent 96 molecule configura¬ 
tions of ice Ih using the TTM3-F model. They are gener¬ 
ated with an algorithm that goes over all possible allowed 
proton configurations and produces structures with no 
net dipole moment^^. 


B. Simulation Procedure 

We used siesta code^i^ to perform DFT calculations 
within the generalized gradient approximation (GGA) to 
the exchange and correlation (XG) functional. The cal¬ 
culations use PBE and vdW-DF^®^ functionalsto 
compare non-local van der Waals effects with semi-local 
GGA approximations. These density functionals have 
previously been shown to give good results for volume 
calculations of hexagonal ice Ihi^ 

Full structural relaxations for calculating the Eq{V) 
curve are performed with t^-l-p basis. For these relax¬ 
ations, we have used a real-space mesh cut-off of 500 Ry 
for the integrals, electronic k-grid cut-off of 10 A (corre¬ 
sponding to 38 k-points) for unit cell calculations of ice 
Ih, force tolerance of 0.001 eV/A , and a density matrix 
tolerance of 10“® electrons. Instead of doing a variable 
cell optimization, we calculate the energy of the relaxed 
structure at a fixed volume for each lattice parameter. 

Even though results with the t^-l-p basis are accurate 
enough to obtain general structural properties'^, for pre¬ 
cise order-disorder free energy values, the energy must be 
very well converged. Recently, a systematic method to 
obtain the finite-range atomic basis sets for liquid water 
and ice has been proposed^. We use the quadruple-^ 
double polarized (qC-l-dp) basis obtained with the new 
proposed framework. We calculate the energy of the 
structures again, with the qC-bdp basis, without relax¬ 
ation. 

For the t^-l-p basis, the error compared to the qC+dp 
basis is —0.23% in lattice constant a, —0.28% in c, and 
—0.71% in the total volume. The change in the energy, 
Eo{Vo) from t^-l-p basis to the qC+dp basis without re¬ 
laxation is 948.6 meV. Further relaxing the structures 
with the qC-l-dp basis does not change the lattice param¬ 
eters, and changes the energy only by 1.3 meV. Details of 
this calculation and the lattice parameters with q^-|-dp 
basis are given in the Supplemental Material (SM)^ . 

For the free energy calculations which include the nu¬ 
clear quantum effects, the vibrational modes are calcu¬ 
lated using the frozen phonon approximation. All the 
force constant calculations are performed with the t^-l-p 
basis. There are two reasons for this: the tC-l-p basis gives 
a good first approximation to the configurational infor¬ 
mation, and the q^-l-dp basis is costly in computer time. 
In addition, the largest error in the free energy calcula¬ 
tions comes from the initial Eq^V) contribution, which we 
reduce significantly as explained above. The error in the 
zero point energy contribution is much smaller than the 
electronic energy error. The force constant calculations 
of proton-ordered ice XI structure use a finer real-space 


mesh cut-off of 800 Ry and an atomic displacement of 
Aa; = 0.06 A. Similarly, the force constant calculations 
of proton-disordered ice Ih use a real-space mesh cut¬ 
off of 500 Ry and an atomic displacement of Ax = 0.08 
A for the frozen phonon calculations. The acoustic sum 
rule has been used throughout the study. 

The phonon frequencies, uJk{Vo) and Griineisen param¬ 
eters 7 fc(Vb) are obtained by diagonalizing the dynamical 
matrix, computed by finite differences from the atomic 
forces in a (3 X 3 X 3) supercell, at volumes slightly be¬ 
low and above Vq. We tested these parameters to ob¬ 
tain force constants in phonon calculations, so that the 
Griineisen parameter calculations have minimum noise^^. 
The Griineisen parameters are calculated for 3 volumes 
corresponding to isotropic expansion and compression 
around the minimum. In order to cover the full Brillouin 
zone of ice XI and ice aXI, 729 k-points are selected, di¬ 
viding each reciprocal lattice vector into 9 equal sections. 


IV. RESULTS 

A. Order-Disorder Phase Transition 

To understand the phase transition between proton- 
disordered and proton-ordered ice, we calculated the co¬ 
hesive energy from eqs. [SI 13 The cohesive energy in¬ 
cluding the zero-point nuclear quantum effects are also 
presented in Table H] The results from different proton- 
disordered configurations using the TTM3-F model all lie 
within ±0.22 meV of each other, as indicated in the first 
line of Table H) The change in the cohesive energy due 
to the residual entropy of hydrogen disorder is on the 
order of 0.22 meV, which means that our quantitative 
prediction of the most stable phase is within this range. 

Both DFT functionals predict stability to decrease in 
the order ice XI —ice aXI —ice Ih. This agrees with 
the experiments that the structure of the ordered phase 
is ferroelectric. On the other hand, TTM3-F predicts 
the stability order to be the reverse, ice Ih —>■ ice aXI 
—^ ice XI, giving the wrong ground state and no phase 
transition. Furthermore, considering the error in the co¬ 
hesive energy, it is impossible to predict the correct sta¬ 
ble phase at the zero temperature limit with this model. 
The phase transition can only be obtained with the DFT 
calculations. 

In order to analyze the proton order-to-disorder phase 
transition temperature, we study the Helmholtz Free en¬ 
ergy at zero pressure. We evaluate the volume depen¬ 
dence of free energy, F(V) at fixed temperature and find 
the value of free energy minimum, E{VF^i^{T)). There¬ 
fore, we obtain a temperature dependence of free energy, 
by evaluating free energy minimum for each temperature, 
using eq. |5]or|4l 

In the classical limit of the free energy, without con¬ 
sidering nuclear quantum effects, as given in eq. [3 DFT 
predicts a phase transition, regardless of the chosen func¬ 
tional. In addition, both semi-local PBE and non-local 
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TABLE I. Classical (E^) and quantum {Ec ) cohesive energies, in meV. Quantum values include zero point effects. 


FF/XC 

Ice 

Eg 

H 2 O 

D 2 O 

H2i®0 

TTM3-F 

Ih 

601.07±0.22 

521.17±0.22 

536.24±0.22 

522.87±0.23 

TTM3-F 

aXI 

600.30 

520.33 

535.41 

522.04 

TTM3-F 

XI 

599.71 

520.11 

535.11 

521.81 

PBE 

Ih 

620.44 

502.07 

526.70 

504.15 

PBE 

aXI 

626.21 

507.17 

531.95 

509.25 

PBE 

XI 

629.06 

509.02 

534.07 

511.11 

vdW-DpP®® 

Ih 

723.94 

601.64 

627.65 

603.64 

vdW-DpP®® 

aXI 

725.53 

602.58 

628.80 

604.57 

vdW-DpP®® 

XI 

728.75 

605.18 

631.52 

607.18 


vdW-DF^®^ functionals overestimate the phase transi¬ 
tion temperature, when nuclear quantum effects are not 
included in the calculations. On the other hand, the 
TTM3-F force field model does not correctly predict the 
stable phase in the low temperature limit, and the differ¬ 
ence between the free energies of the two phases increases 
with temperature. Therefore, it does not show a phase 
transition. We present the full temperature dependence 
of classical free energy in the SM^. To understand how 
each component of eq. [S] contributes to the total free 
energy, we also present the temperature dependence of 
Eo{VF^i^{T)), and —TS terms separately in the SM^. 


B. Isotope Effects in the Transition Temperature 

Going further, the nuclear zero point effects are cal¬ 
culated to compare the predicted phase transition tem¬ 
perature to the experiments. Fig. |4] shows the tempera¬ 
ture dependence of the free energy with zero point effects 
for H 2 O. At all temperatures, TTM3-F model predicts 
ice Ih to be the stable phase. DFT correctly predicts 
the most stable phase as the ferroelectric-ordered ice XI 
for low temperatures, with an energy difference of ^ 4.5 
meV. As the temperature increases, there is a crossing 
at r = 91 K and the proton-disordered ice Ih becomes 
the stable phase beyond this temperature for H 2 O. DFT 
predicts antiferroelectric ice aXI to have lower free en¬ 
ergy than ice Ih at low T, but at all T, ferroelectric 
ice XI is preferred to ice aXL The crossover from the 
antiferroelectric-ordered ice aXI to proton-disordered ice 
Ih is at much lower temperatures, because ice aXI is less 
cohesive than the ferroelectric-ordered ice XL Therefore, 
we establish that with DFT, the most stable phase at low 
temperatures is the ferroelectric-ordered ice XI, in agree¬ 
ment with the experiments. For the rest of the transition 
temperature discussion, we will focus on the ferroelectric- 
ordered to disordered transition, ice Xl/ice Ih. 

Inclusion of zero point effects also allows us to obtain 
the isotope effect in the phase transition temperature, 
since it is experimentally known that the order-disorder 
transition temperature of heavy ice (D 2 O) is higher than 



FIG. 4. Relative free energy per molecule including the quan¬ 
tum zero point effects as a function of temperature. The 
lines show DFT results with vdW-DF^®® functional and the 
dashed lines are the results with the TTM3-F model. DFT 
correctly predicts the most stable phase as ice XI for low tem¬ 
peratures, with an energy difference of ~ 4.5 meV. For low 
temperatures, the TTM3-F model predicts ice Ih as the stable 
phase with ~ 1 meV energy difference at zero temperature; 
and a separation of energy at higher temperatures, making 
the prediction correct for high temperatures only. The re¬ 
sults of the TTM3-F model for ice XI and ice aXI represented 
by black and red dashed lines respectively are almost indis¬ 
tinguishable in this scale. 

light ice (H 2 O) by 4 K— 1 ^. Table ITTl shows that we al¬ 
ready observe the phase transition with calculations at 
the classical limit of free energy, for both PBE and vdW- 
DFPBE approximations, and that the transition temper¬ 
ature decreases with the inclusion of zero-point effects. 
The vdW-DFPBE results are below the glass transition 
temperature around 100-110 K— where proton mobility 
diminishes. This is in general agreement with the ex¬ 
perimental order-disorder phase transition temperatures. 
Although PBE gives a correct prediction of the stable 
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TABLE II. The classical (T°) and quantum (Tc) proton order to disorder transition temperature, Tc (K) including zero point 
effects for ice Ih-ice XI and ice Ih-ice aXI. The ratio of the temperature for different isotopes is given as R(D)=D20/H20 and 
R(^® 0 )=H 2 ^® 0 /H 20 , and the isotope effect on the temperature with respect to the H 2 O transition temperature is also given 
as the isotope effect percentage: IS(A-B)=|J^ — 1. 


Ice 

Method 

rpO 

H 2 O D 2 O H 2 

R(D) R( 

18 O) 

IS(D-H) IS(18 

0 

0 

aXI 

PBE 

153 

151 

156 

151 

1.03 

1.00 

-1-3.31% 

0.00% 

aXI vdW-DpPBE 

42 

30 

35 

30 

1.17 

1.00 

-1-16.67% 

0.00% 

XI 

PBE 

221 

202 

215 

203 

1.06 

1.01 

-1-6.44% 

-10.50% 

XI 

vdW-DpPBE 

105 

91 

97 

90 

1.07 

0.99 

-1-6.59% 

-1.10% 

XI 

ExptiiS 


72 

76 


1.06 


-1-5.56% 



phase, and an isotope effect of 6.4%, the value of the 
phase transition temperature is much larger than the ex¬ 
perimental range. In agreement with the experimental 4 
K difference in the phase transition temperature of the 
isotopes, the vdW-DF^®^ functional predicted transition 
temperature of the heavy ice is larger than the light ice 
with a 6 K difference. As a result, with this method, the 
ratio between the phase transition temperatures of heavy 
and light ice is reproduced within 1% of the experimen¬ 
tal value and the isotope effect on the temperature with 
respect to the H 2 O transition temperature is calculated 
to be 6.6%, as compared to the 5.6% of the experimen¬ 
tal isotope effect. Therefore, it is important to note that 
inclusion of non-local van der Waals forces is critical for 
a reasonable prediction of the transition temperature. 

To understand the main reason behind the differ¬ 
ence between the transition temperatures of different 
isotopes, we study the temperature dependence of each 
component of eq. |4] separately, as presented in Fig. [5] 
The electronic energy difference between the two ices, 
Ao(Ih) — i?o(XI), is larger for H 2 O than D 2 O. This would 
result in a larger transition temperature for H 2 O than 
D 2 O, contradicting experiments. The last two terms, 
zero-point-free vibrational entropy and energy —TSy + 
- Ezp = and 

configurational entropy —TSh also result in larger energy 
difference for H 2 O than D 2 O. However, the zero-point vi¬ 
brational energy Ezp = J2k^k(yF^,^{T))/2 difference 
between the two ices, Ezpijiy) — Azp(XI), is smaller for 
H 2 O than D 2 O. This is the only term that shifts the 
transition temperature of H 2 O below that of D 2 O. These 
results show that the phase transition occurs at a lower 
T for H 2 O than D 2 O because of the zero-point energies 
of the phonon modes. 

This can also be seen from a simple model. The tran¬ 
sition occurs when the free energies of the two ices are 
equal. For the sake of simplicity, we can set the zero of 
energy at the frozen lattice cohesive energy of ice XI and 
denoting by Ed = Ao(Ih) — Ao(XI) and Sp the energy 
and the residual entropy caused by the disorder of ice Hi. 
The free energies are F{X1) = 0 and A(Ih) = Ed — TSp- 
Therefore, at the zeroth order, in the classical limit, it 
follows that Tc(0) = Ed/Sp- When the vibrations are 
included, the free energies from eq. H] are equal at the 



FIG. 5. Top, free energy difference per molecule between ice 
Ih and ice XI calculated with vdW-DF^®® functional in the 
region of the phase transition. Bottom, contributions to this 
free energy difference by each term in eq. |4l Left frozen lattice 
electronic term. Middle zero-point vibrational energy. Right 
remaining terms. All the energies on the bottom plots have 
been shifted to allow them to be compared in the same energy 
scale. 


transition temperature. 


Te = 


Sp 


Ed 


fcsEfcln 


sinh{hujk (Ih) j2kB'Tc) 
sinh(hujk(Xl) / 2kBTc) 


( 8 ) 


A more detailed discussion with the high and low T limits 
of this transition temperature and the latent heat can be 
found in the SM^. Assuming the shifts are not large, 
the isotope shift in the transition temperature can be 
simplified to 


T^D20 rpli20 

^ C ^ C 


^H20 


i?fc(Ih) 


( 9 ) 
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where 


^ ^ sinh(fe^P-0(Ih)/2fcBTj 

sinh(fi4^°^°(Ih)/2fcBTe)’ 


( 10 ) 


and similarly for Rk (XI). Then the low temperature limit 
becomes, 


ATc(low) 

^H20 


h 


^[(cu°=0(Ih)-a;°=0(XI)) 


-(a;«=°(Ih)-.;H^O(XI))] (11) 


where the difference in the frequencies exactly corre¬ 
sponds to the energy difference shown in the bottom mid¬ 
dle panel of Fig. [S] This model clearly shows that the 
main source of the isotope effect in the transition tem¬ 
perature is the difference in the zero point energies of the 
different ices. 



FIG. 6. Vibrational density of states for H 2 O for proton- 
ordered ice XI and disordered Ih structures, as obtained with 
vdW-DF^®® functional. Average Griineisen constants of the 
different modes are given in color code. The inset above zooms 
into the stretching modes and shows the redshift, while the 
inset below zooms into the librational modes and shows the 
blueshift in ice XI with respect to ice Ih. 

This is also evident from the phonon density of states 
of the two ices. Fig. [6] shows phonon density of states 
for H 2 O for both proton-ordered ice XI and proton- 
disordered ice Ih at zero temperature. The colors repre¬ 
sent the average Griineisen parameter of each band sep¬ 
arately. The main effect driving the isotope differences is 
associated to the blue shift of the librational band in ice 
XI with respect to ice Ih, and a corresponding redshift 
of the stretching band. Therefore, with proton ordering, 
the covalency of the intra-molecular bonds is weakened, 
while the inter-molecular hydrogen bonding is strength¬ 
ened. This combined with the weights of the Griineisen 
parameters results in an overall slightly larger zero point 
energy for ice XI than ice Ih. 

One of the reasons of the quantitative difference from 
the experimental results of transition temperature can 
be due to the error in the estimation of residual entropy 
from disorder in both systems. Experimentally, it has 
been shown that at the transition, ice XI loses much but 


not all of the entropy at T^. However, it is not clear 
whether this arises from equilibrium thermal disorder in 
ice XI, or from failure to complete the phase transition, 
leaving some domains of non-equilibrium ice Ih coexist¬ 
ing with ice Xl2^. Another reason for the quantitative 
difference can be the loss of precision of the QHA at 
larger temperatures, as the temperature dependence of 
the phonon vibrations is not taken into account. This 
is also the case in the calculated Vq with isotope effects; 
the calculated values deviate from the experimental val¬ 
ues at larger temperatures^. However, this deviation is 
not significant at around 100 K, which is the region of 
interest of this work. 

Finally, the exact values depend on the choice of the 
DFT functional. While we have shown that the inclu¬ 
sion of vdW interaction in the functional is crucial, it 
should be noted that the local part of the XC functional 
also changes the structure significantly. In Ref. |H, it 
has been shown that vdW-DF functional with the local 
XG flavor of revPBE softens the structure such that the 
anomalous isotope effect on the V(T) is not reproduced 
at low temperatures. In addition. Ref. studied the 
phonon dispersion of ice XL While the distribution of 
the modes are almost identical, the values of the stretch¬ 
ing modes are higher and the librational modes are lower 
by ^ 50 cm“^ than those calculated in this work. Fur¬ 
thermore, Ref. m studied the non-local vdW functionals 
with different GGA and hybrid functionals for the local 
XG, and showed that the cohesive energies depend on 
the choice of these functionals. A hybrid functional with 
exact exchange for the local XC, with a vdW functional 
for the non-local correlations could be a good candidate 
to improve on these results. 

All in all, the QHA within DFT with non-local vdW 
forces, predicts a 6 K temperature difference between the 
isotopes, as compared to the experimental 4 K difference. 
This isotope shift is solely due to the nuclear quantum 
effects from the phonon vibrational energy differences, 
and it is predicted without invoking any other effects, 
such as tunneling. 


V. CONCLUSION 

In this study, we did a detailed analysis of the phase 
transition between the ferroelectric vs. antiferroelectic- 
proton-orderered ice XI and disordered ice Ih. Ab ini¬ 
tio DFT is necessary to correctly predict the most sta¬ 
ble phase of ice as the ferroelectric-ordered ice XL The 
TTM3-F force field model needs improvement for the en¬ 
ergy predictions, especially at low temperatures. 

By including nuclear quantum effects to the free en¬ 
ergy, we have predicted the ferroelectric order to disorder 
phase transition for hexagonal ices. The best accuracy 
requires using the vdW-DF^®® functional, with a tran¬ 
sition temperature at about 91 K for H 2 O and 97 K for 
D 2 O. This 6 K temperature difference is mainly due to 
the difference in the zero point energy of ice with dif- 




























ferent isotopes, while entropy related terms contribute 
in the opposite direction. The method is robust to cor¬ 
rectly predict and explain the isotope effect on the or¬ 
der/disorder phase transition of hexagonal ice Ih. 
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